function [S_EGPP,D_EGPP,P_EGPP,info_EGPP] = EGPP(A_EGPP)
//
// Funzione EGPP: calcola la fattorizzazione LR della matrice A_EGPP.
//
// A_EGPP: matrice N x N;
// S_EGPP, D_EGPP: fattorizzazione LR di P_EGPP A_EGPP (S_EGPP triangolare
//                 inferiore con S_EGPP(k,k) = 1, D_EGPP triangolare superiore,
//                 S_EGPP D_EGPP = P_EGPP A_EGPP);
//
// info_EGPP = -1: la matrice assegnata non è quadrata;
//           =  0: il pivot al passo k è nullo;
//           =  1: la procedura ha determinato la fattorizzazione.
//
// Se info_EGPP ~= 1 si pone S_EGPP = D_EGPP = P_EGPP = %nan.
//
info_EGPP = 1; // resta tale solo se non vi sono errori nel corso della
               // procedura
S_EGPP = %nan; D_EGPP = %nan; // restano tali solo se EG(A_EGPP) non è definita
[NR, NC] = size(A_EGPP);
if NR ~= NC then
   info_EGPP = -1; printf("\n *** Errore: matrice non quadrata.\n\n");
else N = NR;
end;
//
if info_EGPP == 1 then
   P_EGPP = eye(N,N);
   for k = 1:N-1, // colonne da 1 a N-1
       [mx,idx] = max(abs(A_EGPP(k:N,k))); idx = idx + (k-1);
       if idx ~= k then
          riga = P_EGPP(k,:); P_EGPP(k,:) = P_EGPP(idx,:); P_EGPP(idx,:) = riga;
          riga = A_EGPP(k,:); A_EGPP(k,:) = A_EGPP(idx,:); A_EGPP(idx,:) = riga;
       end;
       if A_EGPP(k,k) == 0 then
          info_EGPP = 0; printf("\n *** Errore: pivot nullo al passo %d\n\n",k);
          return;
       else for j = k+1:N, // righe da k+1 a N
                A_EGPP(j,k) = A_EGPP(j,k) / A_EGPP(k,k);
                A_EGPP(j,k+1:N) = A_EGPP(j,k+1:N) - A_EGPP(j,k)*A_EGPP(k,k+1:NC);
            end;
       end;
    end;
    D_EGPP = triu(A_EGPP);
    S_EGPP = eye(A_EGPP) + tril(A_EGPP,-1);
end
endfunction
